function [foF1, NmF1] = foF1_NmF1(foE, foF2)
%
% F1-Layer Critical Frequency
%
%DESCRIPTION:
%This function computes the F1-Layer Critical Frequency (foF1, [MHz]) and
%Maximum Density (NmF1).
%
%PROTOTYPE:
% [foF1, NmF1] = foF1_NmF1(foE, foF2)
%
%--------------------------------------------------------------------------
% INPUTS:
%   foE        [1x1]       E-Layer Crit. Freq.       [MHz]
%   foF2       [1x1]       F2-Layer Crit. Freq.      [MHz] (optional)
%--------------------------------------------------------------------------
% OUTPUTS:
%   foF1       [1x1]       F1-Layer Crit. Freq.      [MHz]
%   NmF1       [1x1]       F1-Layer Max. Density     [1e+11 m-3]
%--------------------------------------------------------------------------
%
%NOTES:
% (none)
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% 09/02/2022 - Updated first if and added third elseif following [2]
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
% [2] "Electron Density Models and Data for Transionospheric Radio
%      Propagation" - Report ITU-R P.2297-1 (05/2019)
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%F1-Layer Critical Frequency [MHz]
if  (foE >= 2) %MHz %[2]
    foF1 = 1.4 * foE;
elseif foE <  2 %MHz
    foF1 = 0;
end
if foF1 < 1e-6
    foF1 = 0;
elseif ((1.4*foE) > (0.85*foF2)) %[2]
    foF1 = 0.85*1.4*foE;
end

%Reduction of 15% of foF1 if too close to foF2 (optional step)
if nargin == 2
    if abs(foF1 - foF2) < 1e-6
        foF1 = foF1 * (1-0.15);
    end
end

%F1-Layer Maximum Density [1e11 m-3]
if (foF1 <= 0) && (foE > 2)
    NmF1 = 0.124 * (foE + 0.5)^2;
else
    NmF1 = 0.124 * foF1^2;
end

end